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Recent experiments show that charge-density wave correlations are prevalent in underdoped 
cuprate superconductors. The correlations are short-ranged at weak magnetic fields but their in¬ 
tensity and spatial extent increase rapidly at low temperatures beyond a crossover field. Here we 
consider the possibility of long-range charge-density wave order in a model of a layered system where 
such order competes with superconductivity. We show that in the clean limit, low-temperature long- 
range order is stabilized by arbitrarily weak magnetic fields. This apparent discrepancy with the 
experiments is resolved by the presence of disorder. Like the field, disorder nucleates halos of charge- 
density wave, but unlike the former it also disrupts inter-halo coherence, leading to a correlation 
length that is always finite. Our results are compatible with various experimental trends, including 
the onset of longer range correlations induced by inter-layer coupling above a characteristic field 
scale. 


I. INTRODUCTION 

The pseudogap state of the cuprate high-temperature 
superconductors (HTSCs) harbors various fluctuating 
electronic orders^ In particular, recent nuclear mag¬ 
netic resonance (NMR) 2 .- and x-ray scattering^SS mea¬ 
surements have found evidence of charge-density wave 
(CDW) fluctuations across this family of materials. The 
observed strength of the CDW fluctuations is anti¬ 
correlated with superconductivity (SC) in the sense that 
the intensity of the CDW scattering peak grows as 
the temperature is reduced towards the superconduct¬ 
ing transition temperature, T c , and then decreases or 
saturates upon entering the SC phase. In addition, the 
CDW signal is enhanced when a magnetic field is used 
to quench SC, while the effect of a magnetic field above 
T c is negligible. Finally, optical excitation of apical oxy¬ 
gen vibrations promotes transient superconducting sig¬ 
natures in YBa2Cu3 06 +a; j 21 ' 22 resembling similar results 
in Lai.675Euo.2Si'o.i25Cu04^ where they were conjec¬ 
tured to be a consequence of the melting of charge stripe 
order 

Motivated by these findings, Hayward et al. 25 i 26 
proposed a phenomenological non-linear sigma model 
(NLSM), which formulates the competition between fluc¬ 
tuating SC and CDW order parameters. Similar models 
emerge also from more microscopic considerations . 2 ^ 30 
Using Monte-Carlo simulations, Hayward et al. calcu¬ 
lated the temperature dependence of the x-ray structure 
factor in the absence of a magnetic field and showed that 
it exhibits a maximum slightly above the Berezinskii- 
Kosterlitz-Thouless temperature, Tbkt , of their two- 
dimensional model. The fact that a similar peak appears 
in zero-field x-ray scattering from YBa2Cu30g-|-ar r — and 
La2- a; Sr a; Cu04j^ was taken as an encouraging sign that 
the NLSM is able to capture salient features of the data. 
The situation, however, is more complicated and the 
structure factor of HgBa2Cu04+^ shows no such peak. 


It is therefore interesting to explore the extent to which 
one can reproduce and understand the various trends 
revealed by experiments from the perspective of a sim¬ 
ple model of competing orders. In particular, we would 
like to ask this question with regard to the transition 
from short-range correlations at low magnetic fields to 
longer range order at high fields, as detected by NMR-, 
ultrasound^! and most recently x-ray scattering^ mea¬ 
surements. 

To this end we incorporate into the NLSM of Ref. [ 2 ^ 
three additional ingredients that are important for com¬ 
parison with experiments, namely, inter-layer couplings, 
a magnetic field and random pinning potentials. We an¬ 
alyze their effects on CDW signatures and ordering ten¬ 
dencies via a large -N approximation, previously used by 
us to study the consequences of thermally excited vor¬ 
tices in the NLSM^ Averages over disorder are calcu¬ 
lated with the replica method, and emphasis is put on 
the low-temperature SC phase where the effects of the 
additional factors are significant. We also present com¬ 
plementary results of Monte-Carlo simulations, which we 
use to study the model beyond the limits of our analytical 
approach. 

We show that in a clean system, without a magnetic 
field, the competition with SC establishes a threshold 
inter-layer coupling for the stabilization of long-range 
CDW order. On the other hand, in the presence of a 
magnetic field, any small inter-layer coupling suffices to 
induce long-range order between the CDW regions which 
nucleate around the cores of the Abrikosov vortices. 
These results are also reflected in the low-temperature 
CDW structure factor of weakly coupled layers. While 
it vanishes linearly with decreasing temperature in the 
field-free system, it diverges when a field is present. Both 
behaviors are inconsistent with the x-ray data. 

In contrast, qualitative agreement with the experimen¬ 
tal phenomenology is obtained when the effects of a ran¬ 
dom pinning potential are taken into account. Since 
favorable disorder configurations nucleate CDW regions 
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that survive the competition with SC down to zero tem¬ 
perature, the structure factor attains a non-zero finite 
value in this limit. This value grows with magnetic 
field, which adds vortices as CDW nucleation centers, 
but true long-range phase order between the CDW re¬ 
gions, predicted by mean-field theory, is avoided due to 
the Imry-Ma argument^ Nevertheless, we find that as 
the field is increased through the failed mean-field tran¬ 
sition, the correlation length is significantly enhanced by 
the effects of inter-layer couplings. At higher tempera¬ 
tures the structure factor exhibits a maximum close to T c , 
which is washed away by both magnetic field and stronger 
disorder strength, while at even higher temperatures it 
becomes magnetic field-independent. The crossover field 
changes little until somewhat below T C: where it diverges. 
However, the increase in the correlation length across 
it diminishes with temperature. The reflection of these 
trends in experiments indicates that both order compe¬ 
tition and disorder are crucial elements in understanding 
the physics of underdoped cuprates. 

In the next section we introduce the model and lay 
out our findings while moving up in complexity from the 
single clean layer to a disordered system of coupled layers. 
The technical derivations of the results are relegated to 
the appendices. In the final section we discuss our work 
in view of experiments done in underdoped cuprates. 


II. MODEL AND RESULTS 
A. The clean, single layer NLSM 

We begin with the model considered by Hayward et al. 
, 25 for a real 6-dimensional order parameter, equivalent to 
a complex SC field T = m + in 2 and two complex CDW 
fields, = 7T3+m4 and = n^+inQ. Here, for the sake 
of simplicity, we disregard quartic and anisotropic CDW 
terms, which appear in Ref. I 25 I and follow our previous 
strategy^ of using a saddle-point approximation for the 
CDW fields, which is formally justified when their num¬ 
ber is large. Thus, we analyze a system described by a 
complex SC field {ip,ip*}, and TV — 2 real CDW fields 
{n a }, where a = 1... TV — 2, whose Hamiltonian is 

H 0 [ip,n a \ = J d 2 r ||(V + 2ieA)ip\ 2 

N—2 

+ t A ( Vn «) 2 +gnl] }, (1) 

a=l 

where p s is the stiffness of the SC order, A p s is the cor¬ 
responding quantity for the CDW components, and gp s 
is the energy density penalty for CDW ordering. We as¬ 
sume that some type of order (SC or CDW) is always 
locally present, in the sense of its amplitude, but that 
the different order parameters compete, as expressed by 


N—2 


( 2 ) 


0 —&Ho 


the constraint 

M 2 + n l = l - 

a=l 

The free energy Fq is given by 

e -/3F 0 _ J x>ip*VipVn 5 ^|V’| 2 + ^2, n a ~ 1 
= J Vip*VipVnVd e~ pH ° 

Xe iPPs f d 2 n»( |b| 2 + E CI «o-l) ; (3) 

where /3 = 1/T is the inverse temperature. In the limit 
TV —> 00 we integrate over n a while assuming that the 
Lagrange multiplier field, d, attains its saddle-point con¬ 
figuration, a = ia. Since we are focusing on signatures 
of the CDW deep inside the SC phase, T -C Tbkt , we 
also assume that the SC fields ip,ip* take their saddle- 
point configurations. Within this approximation^, the 
free energy of the clean layer is given by 

TV — 2 r 1 

(3F 0 = —-—TV In -j3p s (—AV 2 + g + a) 

+\fiPs J d2r [l( V + 2*eA)V>| 2 + a ( \ip\ 2 - l)] , 

(4) 

where the fields ip, ip* and er are determined by the cou¬ 
pled saddle-point equations 


^ = \pp. [-(V + 2 «A) J + a]^ = 0 , 


(5) 


and 


6(3F 0 _ N-2 
Sc r(r) 2 


Tr 


(—AV 2 -f- g + cr) 5 t 


+ 2 VP* (M 2 - 1) = 0, 


( 6 ) 


with S r (r',r") = S(r' — r)S(r" — r). 

We first consider the case of zero magnetic field, where 
the SC field assumes a uniform configuration, 1 p(r) = ipo, 
and a = 0. By substituting this solution in equation d6j) 
we find that 


IV'ol 2 = 1 - 


T 


(7) 


Tmf ’ 

where the mean-field transition temperature is given by 


Tmf 


TV- 2. 
A 

TV- 2 
47tA 


Tr 

In 


(—V 2 + g/X) Ur 

( 32AA 
V 9 a 2 J ' 


( 8 ) 


To obtain the last expression we regularized the theory by 
putting it on a square lattice of spacing a, and assumed 
ga 2 /A <C 1. The CDW structure factor, is defined by 


S ( c 0 = j 2 Jd 2 rd 2 r'e iq (r r > (n a (r)n a (r')) 


(9) 
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where q is measured from the ordering wavevector of n a , 
L 2 is the layer’s area, and (• • •) denotes thermal aver¬ 
aging. We will concentrate on the peak value Scdw = 
S(q = 0), deferring the q dependence to Appendix [A] 
Here, the uniformity of ip and a readily leads to the re¬ 
sult 

Scdw = —> ( 10 ) 

gps 

which vanishes as T —> 0. 

The situation changes upon applying a magnetic field, 
B. The solution of the saddle-point equations m is 
expected to take the form of an Abrikosov lattice of vor¬ 
tices, whose density is determined by the magnetic field. 
Far away from the vortex cores, a = 0 and ip = ipo, just 
as in the zero field case. However, close to the center of 
each vortex, ip vanishes linearly with the distance from 
the vortex center, and a becomes negative. As a result, 
there is a trapped CDW mode inside each core, in ad¬ 
dition to a continuum of scattering modes, which exists 
also without a magnetic field. Using a tight-binding ap¬ 
proximation (see Appendix 0 for these trapped modes 
in equation ©, we can estimate their contribution to 
Scdw ■ The result depends on the order of limits. At 
a low but non-zero temperature, and B —> 0 (more pre¬ 
cisely, when t ~ ge - c Wa^o/B ^ gj 1 / p s ^ we g nc j that 


T 

Scdw = - V A± 

gps 



Ps T M F J 


T 



( 11 ) 


where <po = 7r hc/e is the flux quantum. Here, and 
throughout the paper, we denote by Ai,bi , and Ci var¬ 
ious numerical constants. In the other limit, of a finite 
magnetic field and T —> 0 (when t^> gT/p s ), we obtain 

Scdw = A 3 —JL e c At/9) P sB/T-i/T MF ). ( 12 ) 

tp s gcpo 

Therefore, the structure factor diverges at low tempera¬ 
tures in the presence of a magnetic field. 

In order to go beyond our mean-field results, we per¬ 
formed Monte-Carlo (MC) simulations of the NLSM, in¬ 
corporating the effects of a uniform magnetic field, as 
appropriate for an extreme type-II superconductor. We 
used standard Metropolis updating to study systems on a 
square Lx L lattice, with L ranging from 32 to 200 sites, 
and cylindrical boundary conditions. We present here 
MC results for the experimentally relevant case, N = 6, 
and additionally set A = 1 and ga 2 = 0.3. For these 
parameters, Tbkt/Ps = 0.345. To facilitate comparison 
with the results of Ref. [25j, we present in Fig. QJi the 
structure factor S^ x /a 2 = 2 Scdw/o? as a function of 
T in a clean layer with and without a magnetic field. 
For B = 0, S$ x vanishes linearly as T approaches zero 
in agreement with equation m, but diverges at low T 
for finite B > 0, as in equation ED- Our results are 
generally independent of the system size, except for the 
B > 0,T —> 0 limit, where the diverging S<s> x increases 
with L. 




FIG. 1. Onset of long-range CDW order in the clean system, 
(a) The structure factor of a clean 48 x 48 layer with A = 1 
and ga 2 = 0.3, as function of T for three values of B/Bq, 
where Bo = (po/^Tia 2 . (b) The CDW ordering temperature, 
estimated from a 32 x 32 x 8 multi-layer with the same pa¬ 
rameters, at B = 0 and B = 0.05-Bo. The lines depict the 
solution of the interlayer mean-field condition, equation m- 


B. Clean coupled layers 

Next, we would like to ask whether a weak inter¬ 
layer CDW coupling is sufficient to stabilize long-range 
CDW order. First, let us note that the diverging 
SC susceptibility of each layer at Tbkt implies that 
any weak interlayer Josephson coupling of the form 
psJsc f d 2 r bl’ii > i +1 + H.c.], where i is the layer in¬ 

dex, induces long-range SC order. However, for weak J$c 
the SC transition at T c = Tbkt [l + b±/ ln 2 (6 2 g/Jsc)\ > 
(see Appendix iBl) . has only a small effect on \ip\ and thus 
on the amplitude and ordering tendencies of n. Con¬ 
sequently, we concentrate on the following multi-layer 
Hamiltonian 

H — ^ ( Hq \lpi-, U Ps J _L J d T ^ ■ (13) 

i a,i 
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CDWs on different layers are coupled capacitively. The 
small amplitude of the charge modulation associated with 
the CDW, and its d-wave natur e) 19 ’ 35 imply a weak CDW 
coupling with a complicated real-space structure. We de¬ 
fer the study of the consequences of such structure to a 
future publication and instead treat here the simplest 
model interaction, as expressed in Eq. (fT3l) . In the fol¬ 
lowing we choose J± > 0, although a purely repulsive 
interaction corresponds to J± < 0. However, the two 
cases are related by the transformation n a} i —> (—1 ) l n a ^, 
which reverses the sign of J± but leaves Hq unchanged. 
Consequently, our conclusions regarding the presence of 
long-range CDW order hold for both attractive and re¬ 
pulsive interactions, with the only difference being a 
change in the c-axis ordering wave-vector from 0 to ir/a. 

To estimate the effect of J± we use the interlayer mean- 
field approximatio n 36 ! 37 (see Appendix [Bj which in the 
absence of a held yields the following condition for the 
putative CDW ordering transition 

1 = 2p s J±x{TcDW ), (14) 

expressed in terms of the in-plane CDW susceptibility 
x(T). For a clean system x(T) = Scdw{T)/T and 
equation m implies that condition m can be ful¬ 
filled only if J± > g/2. When this happens uniform 
CDW order is established, and the interlayer coupling 
term in equation m leads to the effective modification 
g —> g — 2J±(n a ) 2 /(n% t ). Hence, for J± > g/2 and T — > 0 
the effective g turns negative, SC disappears and the sys¬ 
tem becomes purely CDW ordered. 

In the presence of a weak magnetic held and at low 
temperatures, the right hand side of the mean-held con¬ 
dition m acquires an additional factor which scales as 
B/gcj) o- However, more important for establishing the 
qualitative difference compared to the held-free case is 
the low-T divergence of Scdw, equation (fl2l) . This 
means that even for J± —> 0, long-range order between 
the CDW regions around the vortex cores does set in 
at Tqdw = c 2 tp s /[gln(c 3 t/J±)\, and coexists with long- 
range SC order. 

Fig. [TJj depicts Tcdw obtained from the onset tem¬ 
perature of the order parameter d 2 rn at i) 2 in 

a clean 32 x 32 x 8 layered system, as a function of the 
interlayer CDW coupling Jj_a 2 . As expected from the 
above mean-held considerations, we hnd a transition for 
Jj_ > g/2 and B = 0, and down to the lowest accessible 
values of Jj_ when B > 0. In addition, however, for J± 
slightly below g/2 and B = 0, we observe a transition 
to an ordered phase which vanishes at a lower critical 
temperature. This behavior can be traced to a maxi¬ 
mum, X'max, in x = Scdw/T, which gives two solutions 
to equation (m in the range l/2p s Xmax < J± < g/2. 

C. Disordered single layer 

The fact that the behavior of the clean system, de¬ 
tailed above, is at odds with experiments motivates us 




FIG. 2. The structure factor of a disordered layer, (a) S^ x /a 2 
of a disordered 48 x 48 layer with A = 1, ga 2 = 0.3 and V 2 a 2 = 
0.075. The inset depicts S'^/a 2 as function of B, for T = 
0.075p s . (b) Sg> x /a 2 of the same layer at B = 0, for various 
levels of disorder strength. The vertical lines depict Tbkt of 
each system, as deduced from the calculated superconducting 
phase stiffness. 


to consider the effects of a random potential, which pins 
the CDW. We begin with the Hamiltonian of a single 
disordered layer 

H 0 [i/,n a -,V a \ = H 0 [ip,n a ] - p s / d 2 r'^V a n a , (15) 

a 

where V a are ind ependent ran dom Gaussian fields satis¬ 
fying V a = 0 and V a {r)Vp{r') = V 2 5 a p6(r — r'), with the 
overline signifying disorder averaging. 

Applying the replica method to the N —> oo limi t 38 ’ 39 
we adapt the saddle-point equations m to the weakly 
disordered case, and calculate the structure factor Scdw 7 
averaged over realizations of the pinning field. In the zero 
field case we find (see Appendix [C]) 

_ T V 2 

Scdw = - 1 — 5 " i ( 16 ) 

9Ps 9 Z 
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T/ps B/Bq 


FIG. 3. Crossover to longer-range CDW order in the disordered system, (a) The crossover line from short-range to longer-range 
CDW order, as deduced from the interlayer mean-field approximation for a 64 x 64 x 16 disordered system with A = 1, ga 2 = 0.3, 
V 2 a 2 = 0.075 and J±a 2 = 0.1. The red dots depict the onset of enhanced field and temperature dependence of the in-plane 
correlation length, (b) The structure factor, averaged over 50 disorder realizations, as function of temperature. The inset show 
the CDW susceptibility of a single layer, used to derive the mean-field crossover line, (c) Magnetic field dependence of the 
structure factor, (d) Magnetic field dependence of the in-plane correlation length. 


which decreases linearly to a finite value as the temper¬ 
ature is reduced to zero. Such behavior reflects the fact 
that due to the random field certain regions assume local 
CDW order even at T = 0. When the system is subject to 
a magnetic field, superconductivity is suppressed inside 
the vortex cores, around which CDW halos are formed. 
As a result, a larger fraction of the system’s area supports 
pinned local CDW order, and Scdw increases. As long 
as t <C ma x(gT/p s , \J gV 2 ) we find that (see Appendix 

El) 


ScDW = -b 

9Ps 

+ Ai 


Yl 

g 2 


/ a 2 

1 A 

„ V 2 ] 

1-1 — 


I T-A 4 — 

l \Ps 

1'mF j 

9 \ 


B 


(17) 


A similar functional form characterizes the T —> 0 spa¬ 
tially averaged Edwards-Anderson order parameter 



Random-field models in the N —» oo limit do not exhibit 
a glass transition^, and qea > 0 for all T. 

In Fig. [Hi we present results for S$ x , averaged over 
60 disorder realizations, for a layer with V 2 a 2 = 0.075. 
The inset depicts its low-temperature B dependence. In 
accordance with our analytical result, equation (fI71) . S® x 
assumes a finite value for T, B —> 0, and grows linearly 
with both T and B. The error bars in our MC results, 
which grow with increasing B and decreasing T, reflect 
the low convergence rates and sensitivity to initial con- 
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FIG. 4. The inter-plane correlation length. The magnetic 
field dependence of £ c for a 64 x 64 x 16 disordered system with 
A = 1, ga 2 = 0.3, V 2 a 2 = 0.075 and J±a 2 = 0.1, averaged 
over 50 disorder realizations. 

ditions which arise in this limit. However, our MC sim¬ 
ulations clearly show that iS$ x does not diverge when 
B > 0, even for temperatures below the range presented 
in the figure. In Fig. [2]o we show 5$ x /a 2 of the same sys¬ 
tem at B = 0 for various disorder strengths. The figure 
also depicts the BKT transition temperature for each sys¬ 
tem, as deduced from the calculated renormalized super¬ 
conducting phase stiffness p s (T ) and the BKT criterion 
Ps{Tbkt) = 2Tbkt/tt. Our results clearly show that the 
peak in the structure factor, which occurs slightly above 
Tbkt , disappears with increasing disorder strength. At 
the same time the zero-temperature value of S'$ x / a 2 in¬ 
creases and approaches a limiting value as global super¬ 
conducting order is suppressed by the disorder. 


D. Coupled disordered layers 

Finally, consider the disordered version of the coupled- 
layer Hamiltonian m, where each plane is described by 
H 0 [i/j,n a -,V], equation (fl5l) . While interlayer mean-field 
approximation predicts a CDW ordering transition at 
TcdWi given by condition m, the Imry-Ma argument^ 
precludes long-range CDW order in a disordered system 
below four dimensions. Nevertheless, for weak disorder 
we expect the failed thermodynamic transition to leave 
its mark in the form of a crossover, which signifies the 
onset of enhanced CDW correlations within and between 
the planes. 

To test this prediction, we start by evaluating the 
mean-field transition temperature from condition m- 
We do so by calculating the disordered-averaged CDW 
susceptibility, y, for a layer with A = 1, ga 2 = 0.3 and 
Va 2 = 0.075, which exhibits an x-ray structure factor 



T/Ps 


FIG. 5. Temperature dependence of the correlation lengths. 
£ a b and £ c for a 64 x 64 x 16 disordered system with A = 1, 
ga 2 = 0.3, V 2 a 2 = 0.075 and J±a 2 = 0.1, averaged over 50 
disorder realizations. 


with a similar temperature dependence to the one mea¬ 
sured at low fields in YBa2Cu306+a;^ (see Fig. [21a). The 
results, presented in the inset of Fig. [3 Jd, show that \ is 
approximately constant at low T and grows linearly with 
B from 1 /p s g, in accord with our large- N analysis (see 
Appendix 0. When combined with equation m this 
implies, for J± < p s g/2, a mean-field transition that oc¬ 
curs at a critical field which is constant over a wide tem¬ 
perature range and then increases rapidly. Such behav¬ 
ior, depicted in Fig. [3}a. for a system with Jj_a 2 =0.1, re¬ 
sembles that of the transition line into a long-range CDW 
phase, as deduced from ultrasound measurements <21 
Next, let us inquire what features of the mean-field 
transition survive the effects of fluctuations. Fig. [3]:) de¬ 
picts the temperature dependence of the structure factor 
(now defined by the average of (n a (r)n a (r')} over the 
three-dimensional system). While the qualitative fea¬ 
tures follow the ones displayed by the two-dimensional 
layer, the coupled layers exhibit, at low T, a rapidly in- 
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creasing S$ x beyond a characteristic field scale, (see Fig. 
Efc). Another relevant signature is displayed by the in¬ 
plane correlation length, £ a b, defined by the inverse half¬ 
width at half maximum of S'$ x (q). Its field dependence, 
shown in Fig. 01, exhibits an inflection point across the 
same magnetic field scale. Recent x-ray measurements!! 
have found a similar behavior in YBa 2 Cu 306 +a;at a range 
of magnetic fields that is comparable to ours if one iden¬ 
tifies the short-distance scale, a, of the NLSM with 2-3 
Cu-Cu spacings. Using also the temperature dependence 
of £ a b and the corresponding data for the inter-plane cor¬ 
relation length, £ c , (see Figs. Q] and [5]) we can map out 
the crossover field in the B — T plane, and find that it 
follows the mean-field transition line. Hence, the follow¬ 
ing picture emerges: At weak magnetic fields the CDW 
correlations are short-ranged and largely confined to the 
planes. However, as the field increases at low tempera¬ 
tures, it combines with the inter-plane coupling to induce 
a crossover to more extended three-dimensional correla¬ 
tions. This crossover becomes sharper with diminishing 
disorder. 


III. DISCUSSION 

We have shown that competition between CDW and 
SC orders in the presence of disorder, can account for 
many of the trends observed in x-ray scattering exper¬ 
iments. Specifically, nucleation of CDW at regions of 
strong attractive disorder makes Scdw attain a finite 
value at T = 0, even for B = 0£ i 15 i 16 This value increases 
linearly with low B *! due to pinned CDW around vor¬ 
tex cores, as seen in scanning-tunneling experiments!! - — . 
At higher temperatures and for weak disorder both 
our simulations and experiments on YBa 2 Cu 306 + a — and 
L^-zSr^CuCU— exhibit a maximum in Scdw close to 
T c , which disappears with increasing B , while there is 
practically no B dependence beyond this temperature^ 
We find that the peak is also washed away with increas¬ 
ing disorder, a fact which may explain its absence in 
HgBa 2 Cu 04 +( 5 — . We note, however, that the predicted 
linear low-T dependence for B = 0, (see Fig. 0 , is re¬ 
flected in some!!, but not all x-ray data<! Another dis¬ 
crepancy exists with Ref. where Scdw was found to 
decrease upon increasing the amount of oxygen disorder. 

Our results also demonstrate a crossover to a regime 
with longer-range CDW correlations at high magnetic 
fields, in accord with x-ray measurements— Moreover, 
the temperature dependence of the crossover field fol¬ 
lows closely the one observed in NMR! and ultrasound!! 
measurements. A commonly advocated scenario!! - — for 
explaining recent quantum oscillations experiments*!! - — 
invokes long-range CDW order as a cause for Fermi- 
surface reconstruction. It is tempting to associate the 
above mentioned regime of enhanced CDW correlations 
with this scenario. However, the spatial extent of CDW 
ordering needed to explain the quantum oscillations ex¬ 
periments remains to be resolved. 
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Appendix A: CDW spectrum in the Abrikosov 
vortex lattice state 


Consider the NLSM for a clean layer with B > 0 and 
T < Tbkt , where an Abrikosov lattice of vortices is ex¬ 
pected to develop. In terms of the orthonormal eigen¬ 
functions, <t> s ( r), and eignevalues, e s , of the operator 


L = -AV 2 + g + cr(r), 
the saddle-point equations take the form 

[-(V + 2ieA(r)) 2 + cr(r)] ^(r) = 0, 


and 


E 




Pps 
N - 2 


[i- IV'WI 2 


(Al) 

(A2) 


(A3) 


We have previously derived an effective Ginzburg- 
Landau theory for the NLSM!!, and showed that the vor¬ 
tex core radius scales at low temperatures as ro ~ g ~ 1 / 2 . 
From Eq. (IA2I) it then follows that inside the core 
a ~ — 1/tq ~ — g , while cr = 0 for r > ro, where 
\ij}\ = IV’ol- Consequently, one expects that in the pres¬ 
ence of a vortex, the spectrum of L consists of a contin¬ 
uum of scattering states with e s > g , and a discrete set of 
bound states with e s < g. Our numerical solution of the 
saddle-point equations confirms these expectations and 
shows that for A = 1 there is a single bound state, <*?o( r ), 
with eigenvalue £q C g, which decays at large distances 
as (/? 0 (r) ~ exp(—r/r 0 ). 

In the presence of a dilute Abrikosov lattice of vor¬ 
tices, i.e., one for which the inter-vortex distance, R , 
obeys R ro, the small overlap between bound states 
in neighboring cores leads to the formation of a tight- 
binding band </>o,k(r). For a square vortex lattice its dis¬ 
persion takes the form 


£o(k) = io — 21 [cos (k x R) + cos (k v R )], (A4) 

with \k X:V \ < 7 t/R, £o — £o — Aso, and 

Ae 0 = -J d 2 rAa R=0 (r)ipl(r), (A5) 


t = - J d 2 np 0 (r)Aa R=Rx {r)ip 0 (r - Rx) 

-Ae o y d 2 np 0 {r)ip 0 {r - Rx). (A6) 


Here, 


Acr R (r) = ^2 ov(r-R'), (A7) 

R'/R 




where oy(r) is the configuration assumed by a in the 
pr esence of a single vortex. Consequently, usi ng R = 
\J(j>o/B, both Aeo and t scale as gexp(—ci\/gfo/B), 
where Ci is a constant that depends on A. 

Under the specified conditions the scattering states still 
form a continuum with e s > g. Since </>o,k(r) vanish 
rapidly between vortices it follows from Eq. (IA3I) that 

£ [i - 1*' 2 - »MI 2 ] • 

sGscattering 

(A8) 

where <5i/>(r) is appreciable only within the cores. There¬ 
fore, 

£ = 4^2 [W 2 - « r >l 2 + WMH ■ 

(A9) 

whose integral over r gives 

• < Aio > 

with a constant C. Evaluating the integral and using 
IV’ol 2 = 1 - T/Tmf , gives 


1 


c 

4k 2 p ‘ 


£ 0 - 4f 


32t eX P 


(l _£_ 

\ T Tmf 


: t < T/p s rl 

:t^>T/p s rl 


n tps (t Tmf) 

(All) 

The effective action for the CDW fields, n a , is of the 
form (f3p s / 2) J d 2 r^2 a n a Ln a , with the result that 


S( q) = 


1 1 



2 


(A12) 


Bloch’s theorem implies that 

</>s( r ) = e lk ' r u„, k (r), (A13) 


with k lying within the magnetic Brillouin zone, and 
u n . k(r + R) = u n ,k(r) for any position R of the Ny 
vortices in the lattice. Therefore, 


Q f \ 1 N* 

s(q) = —TT E ■ 


Pps L 2 ^ £n(q) 


[ d 2 re lQ ' r M„, q /(r) 
J U.C. 


(AM) 

where the integration is over a unit cell of the Abrikosov 
lattice, and where we decomposed q = q' + Q into its 
projection q' to the first Brillouin zone and a reciprocal 
lattice vector Q. 

Due to the spatial integration, the scattering states 
contribution to Scdw = S(q = 0) is dominated by the 
lowest lying extended state with e s ~ g, which is the 
descendent of the k = 0 state of the system with B = 0. 
Therefore, its integral satisfies (1 /L 2 ) | J d 2 rcj> s ( r)| = 1 — 
0(rg/i? 2 ), where the correction is due to its deviations 
from uniformity in the vicinity of the cores. States with 


e s > g provide further contributions of order 0{tq/R 2 ). 
For the localized band we have 


uo,k(r) 


1 

7W 


E 


e ik.(R-r)^ o ( r 


R), 


(A15) 


implying that for R 3> tq 



(A16) 

since ( r ) is normalized and appreciable within r < tq. 
Consequently, the contribution of the band of core states 
to Scdw is of order (ro/R) 2 T/p s eo(k = 0). Using Eq. 
(IA11I) and combing the two contributions, we finally ar¬ 
rive at Eqs. m and m- 

When applying the above considerations to a 
triangular Abrikosov lattice, one needs to take 
into account that for this geometry £o(k) = 

—At [cos 2 (fc x i?/2) + cos(k x R/2) cos(y/3k y R/2) — 1/2] 
with BR 2 = 2<po/\/3. However, these changes do not 
affect the functional dependence of Scdw , but only 
the various numerical constants, which appear in the 
solution. 


Appendix B: The interlayer mean-field 
approximation 

Here, we trade the coupled-layer problem with an effec¬ 
tive single-layer Hamiltonian. For the case of disordered 
CDW-coupled layers, the latter takes the form 

-HmF.i: — Ho ba,i] 

-PsJ± d 2 r ,i( r )( n aj ( r ))MFJ) (Bl) 

ol j=i±l 

where (• • • )mf,i denotes averaging with respect to Hmf,i- 
We are interested in the vicinity of the putative CDW 
ordering temperature, Tcdw , where we would like to 
treat the J± term perturbatively. For the clean system, 
this is justified by the smallness of (n a j(r)) M F,j close to 
Tcdw ■ In the disordered case (n a j(r))MFj is random 
and may be large in regions where the pinning poten¬ 
tial is strong enough to overcome the effects of thermal 
fluctuations. Hence, when disorder is present we assume 
that J±a 2 is small and obtain to leading order 

(n))MF,z = (U'CK,z(r))o 

~\~PsJ .l j d 2 r 'EE X^,i(r,r , )(n^j(r , ))MFj, (B2) 

J p j=i ±l 

where 

Xap,i{ r > r ') 

= ^ [{n a ,i{r)np,i(r ')) 0 - (ria,i( r ))o (WM( r ')> 0 ] - ( B3 ) 

is the in-plane CDW response function, and (• • • )o signi¬ 
fies averaging with respect to Hq. 
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Next, we average equation (IB2I) over the disorder real¬ 
izations. Since the pinning potentials on different layers 
are assumed independent, any correlations between Xa/3,i 
and (ri/3 l i±i)MF,j±i are of order J±. Thus, to lowest order 
in Jj_ we find that 


(n CK (r))MF — ^PsJt 


dVx(r,r') (n Q (r'))MF, (B4) 


where we used (n a ,i(r))o = 0 and the independence of 
(«a,i(r))MF,* and x a /9,i(r,r') = 6 a px(r,r') on i. For 
B = 0 or when disorder or temperature are effective in 
destroying the Abrikosov lattice, translational invariance 
leads to equation d, which expresses the condition for 
the onset of CDW order in term of the in-p lane CDW 
susceptibility x(T) = (1/A 2 ) J d 2 rd 2 r'x( r, r'). 

In a clean system subject to a weak magnetic field 
at low temperatures, the important contribution to the 
integral in equation (IB4I) comes from the vicinity of vor¬ 
tex cores. This implies an approximate condition for the 
transition, similar to equation m, but with the right 
hand side multiplied by a factor that scales as the ratio 
between the core area Tq ~ 1 /</.— and the area of the 
magnetic unit cell R 2 ~ (po/B. 

In the presence of interlayer Josephson coupling and in 
the absence of disorder, the Hamiltonian reads 

H = H 0 [ipj, n a p] - psJsc ( d 2 r ^ [ip*ipj +1 + H.c.], 

i i 

(B5) 

where i is the layer index. The interlayer mean-field 
approximation amounts to replacing H by an effective 
single-layer Hamiltonian of the form 


Hmf = Hq — 2p s Jsc 


d 2 r [ip* (ip) mf + H.c.]. 


(B6) 


We are interested in using this approximation to esti¬ 
mate T c in the multi-layer system. To this end, we calcu¬ 
late (ip) mf- Since it is small in the vicinity of T c we may 
carry out the averaging over Hmf perturbatively in the 
Josephson coupling term. As a result, in the absence of a 
magnetic field and using the fact that (ip(r)ip(r'))o = 0, 
we obtain the following condition for the SC transition 


T c = 2p s J SC j d 2 r(ip*( r'Mr)) 0 . (B7) 


For weak Jsc, T c lies close to Tbkt of a single layer, 
where ( ip*(r')ip(r)) 0 ss \ipo(T c )\ 2 (e^ e ^~ 0 ^ T ^)q. Since for 


T > Tbkt the SC phase correlations decay exponentially 
over the BKT correlation length £, we obtain 

T c = ^ Ps J sc \MTc)\ 2 e(T c ). (B8) 

On a square lattice Tbkt ~ 0-9Ps\iPo(Tbkt)\ 2 , thereby 
establishing, for T c ss Tbkt, a relation between \ipo(T c )\ 2 
and Tbkt■ Fin ally, using BKT cr itical behavior of 
i(T c ) = r 0 exp [b^/T B KT/(T c - T B kt)\ , where b is a con¬ 
stant, we arrive at 


T c = Tbkt 


4 b 2 


In 2 (0.9/47T Jscrf) 


(B9) 


Appendix C: The NLSM with a random pinning 
potential 


We next consider the NLSM of a single layer with in¬ 
dependent Gaussian random potentials, V a . The system 
is described by the action 


S = (3H 0 - Pp s f d 2 r yV a (r)n a (r) 

a 

- / d 2 r J a (r)n a (r) 

O' 

- J d 2 rd 2 r' ^ K a p(r, r')n a (r)np(r'), (Cl) 


to which we have introduced sources in order to calculate 
correlation and response functions in terms of the free 
energy, 


e -/3F = J Xf^VipVn 5 ^| 2 + n 2 - l^J e~ s . 


(C2) 


As a result, 


G a p( r,r') = (n a (r)np( r')) = - 


S/3F 


SK a p( r, r') 


• (C3) 


k =o 


and 


T Xa p(r, r') 


(n a (r)np(r')) - (n a (r))(np(r')) 


S 2 p F 

dJ a (r)SJp( r') 


j=o 


(C4) 


The main difficulty is in calculating F, the free energy 
averaged over realizations of disorder. This can be done 
by employing the replica method in which we consider 
m replicas of the original model. Analytically continuing 
m —» 0 we have F = lim m _>.o F(m)/m, where F(m ) is 
defined by 
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e -/3F(m) = j v^ a *V^ a Vn a a VV a 5 ^|V>°| 2 + £«) 2 - 

= J T>il’ a *'Dip a Fn a T>V a 'Da a e~ ^- ,a 2 ^ / d2r F.c Ha / d2r iffa [l' i / , “l 2 +H a ( n a) 2 - 1 ] _ (C5) 


Integrating over U a and analytically continuing to cr° = ia a , we have e P F b m ) = J Vp’ a *Vvp a Vn , / x V<7 a e with 

1 dp s I (Tr <j 2 ^ |K V + ^eAj^r + <r“WT - 1 JJ +; 

lb o; 


<S(to) = 2^ s J d 2r j£ [1C ^ 7 + 2 ieA)ip a \ 2 + <j a (\ip a \ 2 - 1)] + ££< 

- [ rf2r ££ J «( r KU r ) - [ fl2r(i2r 'EE^(r. r 'K(r)n5( r '), 


<U£ a - /3p s U 2 


(C6) 


ab at 


and L a = — AV 2 + g + cr°(r). Integrating over the CDW fields, n“, gives e P F ( m '> = J V'ip a *'Dip a 'Da a e s ( m \ where 
S'(to) is defined by 


S(to) = iTrln(G _1 - 2 K) + lp Ps J d 2 rJ2 [l(V + 2ieA)^“| 2 + a a (\i/j a \ 2 - 1)] 


1 1 ab 


— 2 / d 2 rd 2 r'J 2 J 2 Ja ( r ) |_( G 1 ~ 2A ') J ( r , r V/3(r') 

ab a(3 ° 


(C7) 


with, 

1 \ab 


We will calculate G ab by assuming a replica-symmetric 
solution of the saddle-point equations, i.e., ip a = ip and 
a a = a. Under this assumption the operator L a = L is 
also replica symmetric, and G ab must obey 


(G-^S^r.r') =pp s [5a b L a -Pp s V 2 ]5 a p6(r-r'), (C8) 
and 

K^{r,r') = 6 ab K a p(r,r'). (C9) 

We would like to calculate the integrals over ip a , , ip a * 
and cr a , using a saddle-point approximation, which is jus- Expanding G aa in the eigen-basis of L 
tified by N —>■ 00 for a a , and provided the disorder is 
weak and satisfies V 2 a 2 = O/l/N), by /3 —> 00 for ip a 
and ip a * , see Eq. (IC22I) below. Within this approxima¬ 
tion, f3F{m ) = S(m), where S(m) is evaluated for the we find the solution 


£ /3ps(S ac L - /3p s V 2 )G cb ( r, r') = S ab S(r - r'). (C14) 


G ab (r, r') = £G^(r)^(r'), (C15) 


configurations of if) 0, ,tp a * and cr a which solve the saddle- 
point equations 


Gft = S st 


Sab 


V 2 


S f3F(m ) 1 a 

5a a (r) = V Ps 


£G““(r,r) + (|r| 2 -l) 


/3p s e s £ s (e s - m/3p s V 2 )\ ’ 


(C16) 


= 0, 


and 


6/3F(m) = U ^_ (v + 2 . gA)2 + ^ = Q (cn) 


expressed in terms of the eigenvalues, e s , of L. 

In the absence of a magnetic field, the saddle-point 
(CIO) equation for ip , Eq. (IC11I) . assumes a uniform solution 
ip = ipo with u = 0. For this case the spectrum of L 
is spanned by plane waves s = k, 0k(r) = e lk r /L and 
£k = Xk 2 + g. Thus, we find for the correlation function 


a k 


Sip a *(r) 2 

Since we are interested in S{m) to order 0(J 2 ,K ), these G a p( r, r') = 5 a p lim i££ Gk°7 2 e ik ' (r " r ' } 

saddle-point equations are not affected by the source ' * * 

fields J and K. Furthermore, since G _1 is diagonal in 
a, /3 and symmetric in a, b and r, r', so is G“^ = 5 a pG ab . = S, 

Thus, we find 


a/3 


f d 2 k 

1 U 2 

1 (2n) 2 

j3p s {Xk 2 + g) (A k 2 +g) 2 _ 


1 xe lk ' (r - r '>, (C17) 

G a p(r, r ) — 5 a p Jiin^ — G (r,r), (C12) from which follows the averaged structure factor 


and, similarly, 


-S'(q) = J 2 J d 2 rd 2 r' e- l ^ r ~^G aa ( r,r') 


A Xa/ 3 (r,r') = 5 a p lim — VG Qb (r,r'). (C13) 

m—>0 TO 


T 


U 2 


ab 


p s (Xq 2 +g) (A q 2 +g) 2 


(CIS) 
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Similarly, we find that the response function is given by 

riMi = 4* ^ E E Gf 

ab k 

1 r d 2 k e ik ( r - r ') 

' Pps J (2tt) 2 Xk 2 +g ’ 


= <5 q 


(C19) 


such that the susceptibility is 

X'(q) = j 2 J d 2 rd 2 r'e~ lcl < ' r - r ' ) Xaa(r,r') 


Ps(Xq 2 +g)' 


(C20) 


Note that unlike the clean case, disorder-induced 
correlations b etween neighboring regio ns lea d to 
(n Q (r))(n Q (r')) ^ 0, and therefore to y(q) ^ S(q)/T. 

Using Eqs. (IC15I) and (IC16I) we obtain that in the 
in, — y 0 limit the saddle-point equation for <r, Eq. (1C 101) . 
takes the form 


f d 2 k 

1 V 2 

(M 

1 (27r) 2 

J3p s {\k 2 +g) (A k 2 +g) 2 _ 

N-2 


which, for A ga 2 , gives 


IV'ol 2 


1- 


(N - 2)V 2 
4nXg 


T 

7^0 ‘ 

1 MF 


(C21) 

(C22) 


where T FIF is the value of Tmf in the clean system. We 
therefore find that disorder reduces IV'ol 2 , as well as Tmf- 
Note that the solution, Eq. (IC22I) . exists only for weak 
enough disorder. For stronger disorder the saddle-point 
configuration is ip = 0 and cr > 0, thus indicating the 
need to take into account fluctuations in ip. 

Next, let us include the effects of a magnetic field on 
Scdw and y. Just as for the clean system, we expect 
that the saddle-point equations of the replicated action 
possess a solution in the form of an Abrikosov lattice. 
Hence, we assume that the spectrum of L consists of a 
continuum of scattering states, similar to those of the 
magnetic-field-free system, and a band originating from 
bound states inside vortex cores. The reasoning that 
was used for the derivation of Eq. (IA9D is then appli¬ 
cable here. If, in addition, we assume that t <C £o, we 
can ignore the dispersion of the tight-binding band and 
approximate it by a flat band with eigenvalue Eo- Con¬ 
sequently, the saddle-point equation for a becomes, as 
m —> 0 


E 

R 


_/3psE 0 


+ 


V 2 


-o J 


l^o(r-R)| 2 

_ [\M 2 - IV>( r )l 2 + I^KOI 2 

~ N - 2 


-. (C23) 


Integrating over r and dividing by the system area gives 


+ 


V 2 


R 2 \_PpsEq e\ 


o J 


C|^o| 2 (^) 2 , (C24) 


where C is a numerical constant. From Eq. (IC24I) we 
find that the assumption t £o is indeed satisfied under 
reasonable conditions, i.e., t <C rna x(gT/p s , gV 2 ). 


In the presence of disorder the expansion of the struc¬ 
ture factor in terms of the eigenstates and eigenvalues of 
L takes the form 




r 1 ^ 2 i 

r 

j3p s Es + E 2 _ 

J d 


(C25) 

Due to the same consideration used for the clean systems, 
we find that the main contribution of the scattering states 
comes from the lowest lying state with e s — g. An ad¬ 
ditional contribution comes from the states cpo( r — R) 
bound to the vortex cores at positions R. Noting that 
L~ 2 J2r 1/ d 2 rip 0 {r — R )| 2 ~ 1 /(gR 2 ), we obtain Eq. 

(ED- 


For the response function one finds 



(C26) 


yielding for T — > 0 and t -C \/gV 2 a q = 0 susceptibility 


1 

p s g 



(po VW 2 ) 


(C27) 


with C a constant. Finally, the T —> 0 spatially aver¬ 
aged Edwards-Anderson order parameter can be easily 
calculated from the saddle-point equation, Eq. (IC23I) . 


Qea(T —> 0) 


J d 2 rG aa { r,r) 

i -\^\ 2 + c \ m 2 (^) 2 


N —2 


(C28) 


which gives Eq. m- 
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